/*! The Matrix library main header file **/
 /** The Matrix library provides a generic data structure to user. Currently this library only supports
  1D, 2D & 3D Matrices.
  The library also provides user with interface for basic Matrix-Scalar, Vector-Vector, Matrix-Vector & Matrix-Matrix
  operations.
  In brief Matrix library allows user for following type of constructors & operations :
  		--CONSTRUCTOR--
  		1) Default constructor.
  		2) Copy Constructor.
  		3) Copy constructor with Arrays.
  		4) Copy Constructor with Initializer lists.
  
  		--OPERATIONS--
  		1) Relational operations.
  		2) Arithmetic operations.
  		3) Scalar operations.
  		4) Matrix-Vector operations.
  		5) Matrix-Matrix operations.
  		6) Copy Assignments.
  		7) Move Assignments.
  		
  		Matrix*Matrix operations is provided only for 1D & 2D matrices.
  
  		Note: We will like to thanks Dr. Bjarne Stroustrup, as this library uses Matrix.h file provided by him.
**/

#ifndef MATRIX_LIB
#define MATRIX_LIB

#include <cstring>
#include <algorithm>
#include <vector>
#include <cstdlib>
#include <time.h>
#include <exception>
#include "Matrix_exception.h"
#include "Matrix_func_obj.h"

using namespace std;

namespace Matrix_lib {
//-----------------------------------------------------------------------------

/*! Index to matrix */
typedef long Index;    

//-----------------------------------------------------------------------------

/*! The general Matrix template is simply a prop for its specializations: */
template<class T = double, int D = 1> class Matrix {
      /*! multidimensional matrix class 	*/
      /*! ( ) does multidimensional subscripting	*/
      /*! [ ] does C style "slicing": gives an N-1 dimensional matrix from an N dimensional one		*/
      /*! row() is equivalent to [ ]*/
      /*! column() is not (yet) implemented because it requires strides.	*/
      /*! slice() to give sub-ranges 	*/
private:
	/*! this should never be compiled	*/
    Matrix();    
};

//-----------------------------------------------------------------------------

template<class T = double, int D = 1> class Row ;    // forward declaration

//-----------------------------------------------------------------------------
//-----------------------------------------------------------------------------

/*! Matrix_base represents the common part of the Matrix classes: */
    /*! matrixs store their memory (elements) in Matrix_base and have copy semantics 	*/
    /*! Matrix_base does element-wise operations	*/
template<class T> class Matrix_base {
public:
	typedef typename vector<T>::iterator iterator;
	/** Default Constructor. **/
    Matrix_base() : sz(0) {}	 

	/** 	Matrix of n elements (default initialized) -- default value of each element set to 0.	**/
    Matrix_base(const Index& n) :sz(n), owns(true), tfer(false)
    {
		elem.assign(n,0);
		iter=elem.begin();
    }
    
    /** 	Matrix of n elements (default initialized) -- default value of each element set to value.	**/
    Matrix_base(const Index& n, const T& value) :sz(n), owns(true), tfer(false)
    { 
		elem.assign(n,value);
		iter=elem.begin();
	}
	
	/** Copy constructor -- matrix of n elements. Every element of type T from Matrix a copied to new matrix object. **/
    Matrix_base(const Matrix_base<T>& p) :elem(p.elem), sz(p.size()), owns(true), tfer(false)
    { iter=elem.begin(); }

    /*! Copy constructor -- matrix of n elements, element from initializer list copied to matrix.  */
    /** For example, if initializer list = {1,2,3,4,5,6}, then matrix has elements (1,2,3,4,5,6) in same sequence. **/
    Matrix_base(const std::initializer_list<T>& a) : elem(a), sz(a.size()), owns(true), tfer(false)
    { iter = elem.begin(); }
    
	/*! Copy constructor -- matrix of n elements, element from all initializer lists inside the initializer list copied to matrix.  */
    /** For example, if initializer list = {{1,2,3} ,{4,5,6}, {7,8,9}}, then matrix has elements (1,2,3,4,5,6,7,8,9) in same sequence. **/
    Matrix_base(const std::initializer_list<std::initializer_list<T>>& a) : sz(a.size()*(*a.begin()).size()), owns(true), tfer(false)
    {
		for(auto x: a)	elem.insert(elem.end(),x.begin(), x.end());		
		iter = elem.begin();
    }

	/*! Copy constructor -- matrix of n elements, element from all initializer lists inside the initializer lists inside initializer list copied to matrix :).  */
    /** For example, if initializer list = {{{1,2,3} ,{4,5,6}}, {{7,8,9}, {10,11,12}}}, then matrix has elements (1,2,3,4,5,6,7,8,9,10,11,12) in same sequence. **/
    Matrix_base(const std::initializer_list<std::initializer_list<std::initializer_list<T>>>& a) : sz(a.size()*(*a.begin()).size()*((*(*a.begin()).begin()).size())), owns(true), tfer(false)
    {
		for(auto x: a)
			for(auto y: x)	elem.insert(elem.end(), y.begin(), y.end());
		iter = elem.begin();
    }
    
	/** Default destructor. **/
    ~Matrix_base()
    { }

	/** Return iterator for first element of Matrix. **/
          iterator begin()       { if (sz==0) throw Matrix_exception("empty vector"); return elem.begin(); }
	/** Return const iterator for first element of Matrix. **/
    const iterator begin() const { if (sz==0) throw Matrix_exception("empty vector"); return elem.begin(); }

	/** Return const iterator for one past last element of Matrix. **/
		iterator end()       { if (sz==0) throw Matrix_exception("empty vector");return elem.end(); }
	/** Return const iterator for one past last element of Matrix. **/
	const iterator end() const { if (sz==0) throw Matrix_exception("empty vector"); return elem.end(); }


	/** Returns size of matrix.	**/
    Index    size() const { return sz; }

	/** Copies individual elements from Matrix a to calling matrix object. **/
    void copy_elements(const Matrix_base& a)
    {
        if (sz!=a.sz) throw Matrix_exception("size not equal");
        for (int i=0; i<size(); ++i) iter[i] = a.elem[i];
    }
    
    /** Copies individual elements from 1D Array of type T a to calling matrix object. **/
    template<int n>
    void copy_elements(const T (&a)[n])
    {
        if (sz!=a.sz) throw Matrix_exception("size not equal");
        for (int i=0; i<size(); ++i) iter[i] = a[i];
    }

	/** Copies individual elements from 2D Array to calling matrix object. **/
    template<int n1, int n2>
    void copy_elements(const T (&a)[n1][n2])
    {
       if (sz!=a.sz) throw Matrix_exception("size not equal");
       for(Index i=0; i<n1; ++i)
		for(Index j=0; j< n2; ++j)
				iter[i*n2+j] = a[i][j];
    }
	
	/** Copies individual elements from 3D Array to calling matrix object. **/
    template<int n1, int n2, int n3>
    void copy_elements(const T (&a)[n1][n2][n3])
    {
       if (sz!=a.sz) throw Matrix_exception("size not equal");
       for(Index i=0; i<n1; ++i)
		for(Index j=0; j< n2; ++j)
			for(Index k=0; k<n3; ++k)
				iter[i*n2*n3+j*n3+k] = a[i][j][k];
    }
    
    /** Interface to copy function to copy elements from Matrix to Matrix **/
    void base_assign(const Matrix_base& a) { copy_elements(a); }

	/** Interface to copy function to copy elements from 1D array to Matrix **/
    template<int n>
    void base_assign(const T (&a)[n]) { copy_elements(a); }

	/** Interface to copy function to copy elements from 2D array to Matrix **/
    template<int n1, int n2>
    void base_assign(const T (&a)[n1][n2]) { copy_elements(a); }

	/** Interface to copy function to copy elements from 3D array to Matrix **/
    template<int n1, int n2, int n3>
    void base_assign(const T (&a)[n1][n2][n3]) { copy_elements(a); }

    /** to get the elements of a local matrix out of a function without copying: (depricated)**/
    void base_xfer(Matrix_base& x)
    {
        if (owns==false) throw Matrix_exception("cannot xfer() non-owner");
        owns = false;     // now the elements are safe from deletion by original owner
        x.tfer = true;    // target asserts temporary ownership
        x.owns = true;
    }
	
	/** Apply function f on elements of Matrix **/
    template<class F> void base_apply(F f) { for (Index i = 0; i<size(); ++i) f(iter[i]); }
    
    /** Apply function f on elements of matrix. Function f takes 1 extra argument to modify value of matrix elements **/
    template<class F> void base_apply(F f, const T& c) { for (Index i = 0; i<size(); ++i) f(iter[i],c); }
    
    
protected:
    vector<T> elem;	/*!< vector is the container which is used to store elements of Matrix. */
    iterator iter; /*!< iterator of type vectorIter (vectorIter is typename for iterator of vector). */
    Index sz;	/*!< Size of Vector container. */
    mutable bool owns;	/*!< This indicates whether the Matrix object owns the storage or not. */
    mutable bool tfer;	/*!< This indicates whether the Matrix object can transfer ownership to other object or not. */
};

//-----------------------------------------------------------------------------

template<class T> class Matrix<T,1> : public Matrix_base<T> {
    Index d1; /*!< dimension of Matrix<T,1>. */

protected:
    /** Used by Row<T,1> */
    Matrix(const Index& n1, Matrix& p) : Matrix_base<T>(n1,p.elem), d1(n1) {}

public:
    /** Default Constructor **/
    Matrix() {}
    
    /** Matrix<T,1> of n elements (default initialized) -- default value of each element set to 0.**/
    Matrix(const Index& n1) : Matrix_base<T>(n1), d1(n1) {}
    
     /** Matrix<T,1> of n elements (default initialized) -- default value of each element set to value.**/
    Matrix(const Index& n1, const T& value) : Matrix_base<T>(n1, value), d1(n1) {}
    
     /** Copy constructor -- Matrix<T,1> of n elements. Every element of type T from 1D Row a copied to new Matrix<T,1> object. **/
    Matrix(const Row<T,1>& a) : Matrix_base<T>(a), d1(a.dim1()) {}
    
    /** Copy constructor -- Matrix<T,1> of n elements. Every element of type T from Matrix<T,1> a copied to new Matrix<T,1> object. **/
    Matrix(const Matrix<T,1>& a) : Matrix_base<T>(a), d1(a.dim1()) { }
    
     /*! Copy constructor -- Matrix<T,1> of n elements, element from initializer list copied to Matrix<T,1>.  */
    /** For example, if initializer list = {1,2,3,4,5,6}, then matrix has elements (1,2,3,4,5,6) in same sequence. **/
    Matrix(const std::initializer_list<T>& a) : Matrix_base<T>(a), d1(a.size()) {}

	/*!Matrix<T,1> of n elements (default initialized) -- elements of 1D array copied to Matrix<T,1>. **/
    template<int n> 
    Matrix(const T (&a)[n]) : Matrix_base<T>(n), d1(n)  { for (Index i = 0; i<n; ++i) this->iter[i]=a[i]; }

	/*! Matrix<T,1> of n elements (default initialized) -- elements of 1D array pointer copied to Matrix<T,1>. */
	/** Note : A user must use a pointer of type T and not an array of type T. Compiler will throw warning for ambigous call if done that. **/
    Matrix(const T* p, const Index& n) : Matrix_base<T>(n), d1(n)  { for (Index i = 0; i<n; ++i) this->iter[i]=p[i]; }

	/*! Construct a new Matrix with element's that are functions of Matrix a's elements: does not modify a unless f has been specifically programmed to modify its argument. */
	/** T f(const T&) would be a typical type for f **/
    template<class F> Matrix(const Matrix& a, F f) : Matrix_base<T>(a.size()), d1(a.d1) { for (Index i = 0; i<this->size(); ++i) this->iter[i] = f(a.elem[i]); }

	/*! Construct a new Matrix with element's that are functions of a's elements: does not modify a unless f has been specifically programmed to modify its argument. **/
	/** T f(const T&, const Arg&) would be a typical type for f **/
    template<class F, class Arg> Matrix(const Matrix& a, F f, const Arg& t1) : Matrix_base<T>(a.size()), d1(a.d1)  { for (Index i = 0; i<this->sz; ++i) this->iter[i] = f(a.elem[i],t1); }
	 
    /** Move constructor -- Matrix now points to same memory as referenced by Matrix a. **/
    Matrix(Matrix<T,1>&& a) {
		xfer(a);
		this->elem=std::move(a.elem);
		this->iter=this->begin();
	}
    /** Move constructor -- Matrix now points to same memory as referenced by Row a. **/
    Matrix(Row<T,1>&& a) {
		xfer(a);
		this->elem=std::move(a.elem);
		this->iter=this->begin();
	}
	
	/** Copy assignment: copy Matrix a into calling Matrix **/
	Matrix& operator=(const Matrix& a);
	
	/** Copy assignment: copy 1D array a into calling Matrix **/
	template<int n>
	Matrix& operator=(const T (&a)[n]);
    
	/** Move assignment : Now both Matrix a and calling matrix share same memory storage **/
	Matrix& operator=(Matrix<T,1>&& a);

	/** Slice operation: Creates a Row with elements from Matrix starting at index n and till m elements from n or end, whichever comes first. **/
	Row<T,1> slice(Index n=0, Index m=0);
	/** Slice operation: Creates a const Row with elements from Matrix starting at index n and till m elements from n or end, whichever comes first. **/
	const Row<T,1> slice(Index n=0, Index m=0) const;

    /** Matrix<T,1> = Matrix<T,1> * Matrix<T,1> **/
    Matrix<T,2> operator*(const Matrix<T,1>& m1);
    /** Matrix<T,1> = Matrix<T,1> + Matrix<T,1> **/
    Matrix operator+(const Matrix<T,1>& m1);
    /** Matrix<T,1> = Matrix<T,1> - Matrix<T,1> **/
    Matrix operator-(const Matrix<T,1>& m1);
    /** Matrix<T,1> = Matrix<T,1> & Matrix<T,1> **/
    Matrix operator&(const Matrix<T,1>& m1);
    /** Matrix<T,1> = Matrix<T,1> | Matrix<T,1> **/
    Matrix operator|(const Matrix<T,1>& m1);
    /** Matrix<T,1> = Matrix<T,1> ^ Matrix<T,1> **/
    Matrix operator^(const Matrix<T,1>& m1);
    
    /** Check if element accessed falls within range of Matrix **/
    void range_check(const Index& n1) const
    { if (n1<0 || d1<=n1) throw Matrix_exception("1D range error: dimension 1"); }

    /** Subscript operation: Returns element at index n, if within range  **/
          T& operator()(Index n)       { range_check(n); return this->iter[n]; }
    const T& operator()(Index n) const { range_check(n); return this->iter[n]; }

	/** Subscript operation: Returns element at index n, if within range  **/
          T& operator[](Index n)       { return row(n); }
    const T& operator[](Index n) const { return row(n); }

	/** Checks if element accessed within range, returns element at index n.  **/
          T& row(Index n)       { range_check(n); return this->iter[n]; }
    const T& row(Index n) const { range_check(n); return this->iter[n]; }
    
	/*! Element-wise assignment operation **/
	/** Note: All other element-wise operations are implemented as a common interface for Matrices of any dimenions **/
	Matrix&  operator=(const T& c)  { this->base_apply(Assign<T>(),c); return *this; }


	/** Transfer Ownership operation (depricated)**/
    void xfer(Matrix<T,1>& a)    // make an Matrix to move elements out of a scope
    {
		d1=a.dim1();
		this->sz=a.size();
		this->owns=a.owns;
		this->tfer=a.tfer;
		a.sz=a.d1=0;
		a.owns=a.tfer=false;
    }
	/** returns dimension of Matrix<T,1> **/
    Index dim1() const { return d1; }
	
	/** Default destructor **/
    ~Matrix() { }
};

//-----------------------------------------------------------------------------

template<class T> class Matrix<T,2> : public Matrix_base<T> {
    Index d1;	/*!< dimension 1 (number of rows) of Matrix<T,2>. */
    Index d2;	/*!< dimension 2 (number of columns) of Matrix<T,2>. */

protected:    
    /** Used by Row<T,2> */
    Matrix(const Index& n1, const Index& n2, const Matrix& p) : Matrix_base<T>(n1*n2,p), d1(n1), d2(n2) {}

public: 
     /** Default Constructor **/
    Matrix() {}

	/** Matrix<T,2> of n1*n2 elements (default initialized) -- default value of each element set to 0.**/
    Matrix(const Index& n1, const Index& n2) : Matrix_base<T>(n1*n2), d1(n1), d2(n2) {}
    
    /** Matrix<T,2> of n1*n2 elements (default initialized) -- default value of each element set to value.	**/
    Matrix(const Index& n1, const Index& n2, const T& value) : Matrix_base<T>(n1*n2, value), d1(n1), d2(n2) {}
  
    /** Copy constructor -- Matrix<T,2> of n elements. Every element of type T from Row<T,2> a copied to new Matrix<T,2> object. **/
    Matrix(const Row<T,2>& a) : Matrix_base<T>(a), d1(a.dim1()), d2(a.dim2()) {}
    
    /** Copy constructor -- Matrix<T,2> of a.dim1()*a.dim2() elements. Every element of type T from 2D Matrix a copied to new Matrix<T,2> object. **/
    Matrix(const Matrix<T,2>& a) : Matrix_base<T>(a), d1(a.dim1()), d2(a.dim2()) {}
    
     /*! Copy constructor -- Matrix<T,2> of n elements, element from initializer list copied to Matrix<T,2>.  */
    /** For example, if initializer list = {1,2,3,4,5,6}, then matrix has elements (1,2,3,4,5,6) in same sequence. **/
    /** d1 = 1, d2 = size of initializer list. **/
    Matrix(const std::initializer_list<T>& a) : Matrix_base<T>(a), d1(1), d2(a.size()) {}

    /*! Copy constructor -- Matrix<T,2> of n elements, element from all initializer lists inside the initializer list copied to Matrix<T,2>.  */
    /** For example, if initializer list = {{1,2,3} ,{4,5,6}, {7,8,9}}, then matrix has elements ({1,2,3} , {4,5,6}, {7,8,9}) in same sequence. **/
    Matrix(const std::initializer_list<std::initializer_list<T>>& a) : Matrix_base<T>(a), d1(a.size()), d2((*a.begin()).size()) {}

	/*!Matrix<T,2> of n elements (default initialized) -- elements of 2D array copied to Matrix<T,2>. **/
    template<int n1, int n2> 
    Matrix(const T (&a)[n1][n2]) : Matrix_base<T>(n1*n2), d1(n1), d2(n2)
    {
        for (Index i = 0; i<n1; ++i)
            for (Index j = 0; j<n2; ++j) this->iter[i*n2+j]=a[i][j];
    }

	/*! Construct a new Matrix<T,2> with element's that are functions of Matrix<T,2> a's elements: does not modify a unless f has been specifically programmed to modify its argument. */
	/** T f(const T&) would be a typical type for f **/
    template<class F> Matrix(const Matrix& a, F f) : Matrix_base<T>(a.size()), d1(a.d1), d2(a.d2)
    {
        for (Index i = 0; i<this->size(); ++i) this->iter[i] = f(a.iter[i]); 
    }

	/*! Construct a new Matrix<T,2> with element's that are functions of Matrix<T,2> a's elements: does not modify a unless f has been specifically programmed to modify its argument. **/
	/** T f(const T&, const Arg&) would be a typical type for f **/
    template<class F, class Arg> Matrix(const Matrix& a, F f, const Arg& t1) : Matrix_base<T>(a.size()), d1(a.d1), d2(a.d2)
    {
        for (Index i = 0; i<this->size(); ++i) this->iter[i] = f(a.iter[i],t1); 
    }
    
    /** Move constructor -- Matrix<T,2> now points to same memory as referenced by Matrix<T,2> a. **/
    Matrix(Matrix<T,2>&& a) {
		xfer(a);
		this->elem=std::move(a.elem);
		this->iter=this->begin();
	}

    /** Move constructor -- Matrix<T,2> now points to same memory as referenced by Row<T,2> a. **/
    Matrix(Row<T,2>&& a) {
		xfer(a);
		this->elem=std::move(a.elem);
		this->iter=this->begin();
	}

	/** Copy assignment: copy Matrix<T,2> a into calling Matrix<T,2> **/
    Matrix& operator=(const Matrix& a);
        
	/** Move assignment : Now both Matrix<T,2> a and calling Matrix<T,2> share same memory storage **/
    Matrix& operator=(Matrix<T,2>&& a);
	
	/** Slice operation: Creates a Row<T,2> with elements from Matrix<T,2> starting at index (n1,n2) and till length1 elements from n1 and length2 elements from n2 or end, whichever comes first. **/
	Row<T,2> slice(Index n1=0, Index n2=0, Index length1=0, Index length2=0);
	
	/** Slice operation: Creates a const Row<T,2> with elements from Matrix<T,2> starting at index (n1,n2) and till length1 elements from n1 and length2 elements from n2 or end, whichever comes first. **/
	const Row<T,2> slice(Index n1=0, Index n2=0, Index length1=0, Index length2=0) const;

    /** Matrix<T,2> = Matrix<T,2> * Matrix<T,2> **/
    Matrix operator*(const Matrix<T,2>& m1);
    /** Matrix<T,2> = Matrix<T,2> + Matrix<T,2> **/
    Matrix operator+(const Matrix<T,2>& m1);
    /** Matrix<T,2> = Matrix<T,2> - Matrix<T,2> **/
    Matrix operator-(const Matrix<T,2>& m1);
    /** Matrix<T,2> = Matrix<T,2> & Matrix<T,2> **/
    Matrix operator&(const Matrix<T,2>& m1);
    /** Matrix<T,2> = Matrix<T,2> | Matrix<T,2> **/
    Matrix operator|(const Matrix<T,2>& m1);
    /** Matrix<T,2> = Matrix<T,2> ^ Matrix<T,2> **/
    Matrix operator^(const Matrix<T,2>& m1);
    
	/** Check if element accessed falls within range of Matrix<T,2> **/
    void range_check(Index n1, Index n2) const
    {
        if (n1<0 || d1<=n1) throw Matrix_exception("2D range error: dimension 1");
        if (n2<0 || d2<=n2) throw Matrix_exception("2D range error: dimension 2");
    }

	/** Slice operation: Returns reference to reference to entire Row<T,1> at row index n1, if within range  **/
          Row<T,1> operator[](Index n1)       { return row(n1); }	//(return a row):

	/** Slice operation: Returns const reference to const reference to entire Row<T,1> at row index n1, if within range  **/
	const Row<T,1> operator[](Index n1) const { return row(n1); } //(return a row):

	/** Checks if element accessed within range, returns reference to entire row at index n1.  **/
          Row<T,1> row(Index n)       { range_check(n,0); return Row<T,1>(n,(*this)); }
          
	/** Checks if element accessed within range, returns rconstant reference entire row at index n1.  **/
    const Row<T,1> row(Index n) const { range_check(n,0); return Row<T,1>(n,(*this)); }


	/** Subscript operation: Returns reference to element at index (n1,n2), if within range  **/
          T& operator()(Index n1, Index n2)       { range_check(n1,n2); return this->iter[n1*d2+n2]; }
	/** Subscript operation: Returns const reference to element at index (n1,n2), if within range  **/
    const T& operator()(Index n1, Index n2) const { range_check(n1,n2); return this->iter[n1*d2+n2]; }

    
	/*! Element-wise assignment operation **/
	/** Note: All other element-wise operations are implemented as a common interface for Matrices of any dimenions **/
	Matrix&  operator=(const T& c)  { this->base_apply(Assign<T>(),c); return *this; }

	/** returns dimension1 of Matrix<T,2> **/
    Index dim1() const { return d1; }
    
    /** returns dimension1 of Matrix<T,2> **/
    Index dim2() const { return d2; }

	/** Transfer Ownership operation (depricated)**/
    void xfer(Matrix<T,2>& a)    // make an Matrix to move elements out of a scope
    {
		d1=a.dim1();
		d2=a.dim2();
		this->sz=a.size();
		this->owns=a.owns;
		this->tfer=a.tfer;
		a.sz=a.d1=a.d2=0;
		a.owns=a.tfer=false;
    }
	
	/** Default destructor **/
    ~Matrix() { }
};

//-----------------------------------------------------------------------------

template<class T> class Matrix<T,3> : public Matrix_base<T> {
    Index d1;	/*!< dimension 1 (number of rows) of Matrix<T,3>. */
    Index d2;	/*!< dimension 2 (number of columns) of Matrix<T,3>. */
    Index d3;	/*!< dimension 3 (number of columns in dimension 2) of Matrix<T,3>. */

protected:
	/** Used by Row<T,3> */
    Matrix(Index n1, Index n2, Index n3, Matrix<T,3>& p) : Matrix_base<T>(n1*n2*n3,p), d1(n1), d2(n2), d3(n3)  {}

public: 
    /** Default Constructor **/
    Matrix() {}
    
    /** Matrix<T,3> of n1*n2*n3 elements (default initialized) -- default value of each element set to 0.**/
    Matrix(Index n1, Index n2, Index n3) : Matrix_base<T>(n1*n2*n3), d1(n1), d2(n2), d3(n3) { }
    
    /** Matrix<T,3> of n1*n2*n3 elements (default initialized) -- default value of each element set to value.	**/
    Matrix(Index n1, Index n2, Index n3, T value) : Matrix_base<T>(n1*n2*n3, value), d1(n1), d2(n2), d3(n3) { }

    /** Copy constructor -- Matrix<T,3> of n elements. Every element of type T from Row<T,3> a copied to new Matrix<T,3> object. **/
	Matrix(const Row<T,3>& a) : Matrix_base<T>(a), d1(a.dim1()), d2(a.dim2()), d3(a.dim3()) {}

	/** Copy constructor -- Matrix<T,3> of a.dim1()*a.dim2()*a.dim3() elements. Every element of type T from Matrix<T,3> a copied to new Matrix<T,3> object. **/
    Matrix(const Matrix<T,3>& a) : Matrix_base<T>(a), d1(a.d1), d2(a.d2), d3(a.d3) {}

	/*! Copy constructor -- Matrix<T,3> of n elements, element from initializer list copied to Matrix<T,3>.  */
	/** For example, if initializer list = {1,2,3,4,5,6}, then matrix has elements (1,2,3,4,5,6) in same sequence with dimension 1 equal to 1. **/
	/** d1 = 1, d2 = 1, d3 = size of initializer list. **/
	Matrix(const std::initializer_list<T>& a) : Matrix_base<T>(a), d1(1), d2(1), d3(a.size()) {}
   
	/*! Copy constructor -- Matrix<T,3> of n elements, element from all initializer lists inside the initializer lists inside initializer list copied to Matrix<T,3> :).  **/
	/** For example, if initializer list = {{{1,2,3} ,{4,5,6}}, {{7,8,9}, {10,11,12}}}, then matrix has elements  accessed as ({{{1,2,3} ,{4,5,6}}, {{7,8,9}, {10,11,12}}}) in same sequence. **/
    Matrix(const std::initializer_list<std::initializer_list<std::initializer_list<T>>>& a) : Matrix_base<T>(a), d1(a.size()), d2((*a.begin()).size()), d3((*(*a.begin()).begin()).size()) {}

	/*!Matrix<T,3> of n1*n2*n3 elements (default initialized) -- elements of 3D array copied to Matrix<T,2>. **/
    template<int n1, int n2, int n3> 
    Matrix(const T (&a)[n1][n2][n3]) : Matrix_base<T>(n1*n2*n3), d1(n1), d2(n2), d3(n3)
        // deduce "n1", "n2", "n3" (and "T"), Matrix_base allocates T[n1*n2*n3]
    {
        for (Index i = 0; i<n1; ++i)
            for (Index j = 0; j<n2; ++j)
                for (Index k = 0; k<n3; ++k)
                    this->iter[i*n2*n3+j*n3+k]=a[i][j][k];
    }

	/*! Construct a new Matrix<T,3> with element's that are functions of Matrix<T,3> a's elements: does not modify a unless f has been specifically programmed to modify its argument. */
	/** T f(const T&) would be a typical type for f **/
    template<class F> Matrix(const Matrix& a, F f) : Matrix_base<T>(a.size()), d1(a.dim1()), d2(a.dim2()), d3(a.dim3())
    {
        for (Index i = 0; i<this->size(); ++i) this->iter[i] = f(a.iter[i]); 
    }

	/*! Construct a new Matrix<T,3> with element's that are functions of Matrix<T,3> a's elements: does not modify a unless f has been specifically programmed to modify its argument. **/
	/** T f(const T&, const Arg&) would be a typical type for f **/
    template<class F, class Arg> Matrix(const Matrix& a, F f, const Arg& t1) : Matrix_base<T>(a.size()), d1(a.dim1()), d2(a.dim2()), d3(a.dim3())
    {
        for (Index i = 0; i<this->size(); ++i) this->iter[i] = f(a.iter[i],t1); 
    }
    
    /** Move constructor -- Matrix<T,3> now points to same memory as referenced by Matrix<T,3> a. **/
    Matrix(Matrix<T,3>&& a) {
		xfer(a);
		this->elem=std::move(a.elem);
		this->iter=this->begin();
	}
    /** Move constructor -- Matrix<T,3> now points to same memory as referenced by Row<T,3> a. **/
    Matrix(Row<T,3>&& a) {
		xfer(a);
		this->elem=std::move(a.elem);
		this->iter=this->begin();
	}


	/** Copy assignment: copy Matrix<T,3> a into calling Matrix<T,2> **/
    Matrix& operator=(const Matrix& a);
        
	/** Move assignment : Now both Matrix<T,3> a and calling Matrix<T,3> share same memory storage **/
    Matrix& operator=(Matrix<T,3>&& a);
    
    /** Slice operation: Creates a Row<T,3> with elements from Matrix<T,3> starting at index (n1,n2,n3) and till length1 elements from n1, length2 elements from n2 and length3 elements from n3 or end, whichever comes first. **/
    Row<T,3> slice(Index n1=0, Index n2=0, Index n3=0, Index length1=0, Index length2=0, Index length3=0);
    
    /** Slice operation: Creates a const Row<T,3> with elements from Matrix<T,3> starting at index (n1,n2,n3) and till length1 elements from n1, length2 elements from n2 and length3 elements from n3 or end, whichever comes first. **/
    const Row<T,3> slice(Index n1=0, Index n2=0, Index n3=0, Index length1=0, Index length2=0, Index length3=0) const;

	/** Slice operation: Returns reference to Row<T,2> at row index n1, if within range  **/
          Row<T,2> operator[](Index n)       { return row(n); }
          
    /** Slice operation: Returns const reference to Row<T,2> at row index n1, if within range  **/
    const Row<T,2> operator[](Index n) const { return row(n); }

	/** Checks if element accessed within range, returns reference to entire Row<T,2> at index n1.  **/
          Row<T,2> row(Index n)       { range_check(n,0,0); return Row<T,2>(n,(*this)); }
    
    /** Checks if element accessed within range, returns const reference to entire Row<T,2> at index n1.  **/
    const Row<T,2> row(Index n) const { range_check(n,0,0); return Row<T,2>(n,(*this)); }

	/** Check if element accessed falls within range of Matrix<T,3> **/
    void range_check(Index n1, Index n2, Index n3) const
    {
        // std::cerr << "range check: (" << d1 << "," << d2 << "): " << n1 << " " << n2 << "\n";
        if (n1<0 || d1<=n1) throw Matrix_exception("3D range, dimension 1");
        if (n2<0 || d2<=n2) throw Matrix_exception("3D range, dimension 2");
        if (n3<0 || d3<=n3) throw Matrix_exception("3D range, dimension 3");
    }

	/** Subscript operation: Returns reference to element at index (n1,n2,n3), if within range  **/
          T& operator()(Index n1, Index n2, Index n3)       { range_check(n1,n2,n3); return this->iter[d2*d3*n1+d3*n2+n3]; }; 
	/** Subscript operation: Returns const reference to element at index (n1,n2,n3), if within range  **/
    const T& operator()(Index n1, Index n2, Index n3) const { range_check(n1,n2,n3); return this->iter[d2*d3*n1+d3*n2+n3]; };

	/*! Element-wise assignment operation **/
	/** Note: All other element-wise operations are implemented as a common interface for Matrices of any dimenions **/
	Matrix&  operator=(const T& c)  { this->base_apply(Assign<T>(),c); return *this; }

    /** Matrix<T,3> = Matrix<T,3> + Matrix<T,3> **/
    Matrix operator+(const Matrix<T,3>& m1);
    /** Matrix<T,3> = Matrix<T,3> - Matrix<T,3> **/
    Matrix operator-(const Matrix<T,3>& m1);
    /** Matrix<T,3> = Matrix<T,3> & Matrix<T,3> **/
    Matrix operator&(const Matrix<T,3>& m1);
    /** Matrix<T,3> = Matrix<T,3> | Matrix<T,3> **/
    Matrix operator|(const Matrix<T,3>& m1);
    /** Matrix<T,3> = Matrix<T,3> ^ Matrix<T,3> **/
    Matrix operator^(const Matrix<T,3>& m1);

	/** returns dimension1 of Matrix<T,3> **/
    Index dim1() const { return d1; }
    
    /** returns dimension1 of Matrix<T,3> **/
    Index dim2() const { return d2; }

    /** returns dimension1 of Matrix<T,3> **/
    Index dim3() const { return d3; }

    // Column<T,2> column(Index n); // not (yet) implemented: requies strides and operations on columns
	/** Transfer Ownership operation (depricated)**/
    void xfer(Matrix<T,3>& a)    // make an Matrix to move elements out of a scope
    {
		d1=a.dim1();
		d2=a.dim2();
		d3=a.dim3();
		this->sz=a.size();
		this->owns=a.owns;
		this->tfer=a.tfer;
		a.sz=a.d1=a.d2=a.d3=0;
		a.owns=a.tfer=false;
    }

    /** Deafult destructor  **/
    ~Matrix() { }
};

//-----------------------------------------------------------------------------

// The default values for T and D have been declared before.
template<class T, int D> class Row {
    // general version exists only to allow specializations
private:
        Row();
};

//-----------------------------------------------------------------------------

template<class T> class Row<T,1> : public Matrix<T,1> {
public:
	Row() : Matrix<T,1> () {}
	
    Row(const Index& n, const Index& length, const Matrix<T,1>& a) : Matrix<T,1>(length)
    {
		for(int i=n,j=0;j<length;++i, ++j) (*this)[j] = a[i];
    }

    Row(const Row<T,1>& a) : Matrix<T,1>(a) {}

    Row(const Matrix<T,1>& a) : Matrix<T,1>(a) {}
    
    Row(const Index& n, const Matrix<T,2>& a) : Matrix<T,1>(a.dim2()) {
		for(Index i=0;i<a.dim2();++i)	(*this)(i)=a(n,i);
	}

    Row<T,1>& operator=(const T& c) {this->base_apply(Assign<T>(),c); return *this; }

    Row<T,1>& operator=(const Matrix<T,1>& a)
    {
		if (this->size()!=a.size()) throw Matrix_exception("length error in Row 1D");
		for(int i=0;i<(*this).size();i++) (*this)[i] = a[i];
        return (*this);
	}
};

//-----------------------------------------------------------------------------

template<class T> class Row<T,2> : public Matrix<T,2> {
public:
	Row() : Matrix<T,2> () {}
	
    Row(Index n1, Index n2, Index length1, Index length2, const Matrix<T,2>& a) : Matrix<T,2>(length1,length2)
    {
		int r=0,s=0;
		for(int i=0;i<length1;++i) 
			for(int k=0;k<length2;++k) {
				r=n1+i;
				s=n2+k;
				(*this)(i,k) = a(r,s);
			}
    }

    Row(const Row<T,2>& a) : Matrix<T,2>(a) {}

    Row(const Matrix<T,2>& a) : Matrix<T,2>(a) {}
    
    Row(const Index& n, const Matrix<T,3>& a) : Matrix<T,2>(a.dim2(),a.dim3()) {
		for(Index i=0;i<a.dim2();++i)
			for(Index k=0;k<a.dim3();++k)
				(*this)(i,k)=a(n,i,k);
	}
	
    Matrix<T,2>& operator=(const T& c) { this->base_apply(Assign<T>(),c); return *this; }

    Matrix<T,2>& operator=(const Matrix<T,2>& a)
    {
		if (this->size()!=a.size()) throw Matrix_exception("length error in Row 2D");
		for(int i=0;i<(*this).dim1();++i) 
			for(int k=0;k<(*this).dim2();++k) {
				(*this)(i,k) = a(i,k);
			}
        return (*this);
    }
};

//-----------------------------------------------------------------------------

template<class T> class Row<T,3> : public Matrix<T,3> {
public:
	Row() : Matrix<T,3> () {}

    Row(Index n1, Index n2, Index n3, Index length1, Index length2, Index length3, const Matrix<T,3>& a) : Matrix<T,3>(length1,length2,length3)
    {
		int q=0,r=0,s=0;
		for(int i=0;i<length1;++i) 
			for(int j=0;j<length2;++j)
				for(int k=0;k<length3;++k) {
					q=n1+i;
					r=n2+j;
					s=n3+k;
					(*this)(i,j,k) = a(q,r,s);
				}
    }

    Row(const Row<T,3>& a) : Matrix<T,3>(a) {}

    Row(const Matrix<T,3>& a) : Matrix<T,3>(a) {}
    
    Matrix<T,3>& operator=(const T& c) { this->base_apply(Assign<T>(),c); return *this; }

    Matrix<T,3>& operator=(const Matrix<T,3>& a)
    {
		if (this->size()!=a.size()) throw Matrix_exception("length error in Row 3D");
		for(int i=0;i<(*this).dim1();++i) 
			for(int j=0;j<(*this).dim2();++j)
				for(int k=0;k<(*this).dim3();++k) {
					(*this)(i,j,k) = a(i,j,k);
			}
        return (*this);
    }
};

/******************************* Matrix<T,1> **********************************/

template<class T>
Matrix<T,1>& Matrix<T,1>::operator=(const Matrix<T,1>& a)
// copy assignment: let the base do the copy
{
	if((*this).size()==0) throw Matrix_exception("copying into uninitialized matrix in Matrix<T,1> operator=");
	if(a.size()==0) throw Matrix_exception("copying from uninitialized matrix in Matrix<T,1> operator=");
	if (d1!=a.d1) throw Matrix_exception("length error in Matrix<T,1> operator=");
	this->base_assign(a);
	return *this;
}

template<class T>
template <int n>
Matrix<T,1>& Matrix<T,1>::operator=(const T (&a)[n])
// copy assignment: let the base do the copy
{
	if((*this).size()==0) throw Matrix_exception("copying into uninitialized matrix in Matrix<T,1> operator=");
	if (d1!=n) throw Matrix_exception("length error in 1D in Matrix<T,1> operator=");
	for(int i=0;i<n;++i) (*this)[i] = a[i];
	return *this;
}

template<class T>
Matrix<T,1>& Matrix<T,1>::operator=(Matrix<T,1>&& a)
// move assignment
{
	xfer(a);
	(*this).elem=std::move(a.elem);
	(*this).iter=(*this).elem.begin();
	return *this;
}

template<class T>
Row<T,1> Matrix<T,1>::slice(Index n, Index m)
// m elements starting with a[n]
{
	if (n<0) n=0;
	else if(d1<n) n=d1;    // one beyond the end
	if (m<=0 || (d1<n+m)) m=d1-n;
	Row<T,1> z(n,m,*this);
	return z;
}

template<class T>
const Row<T,1> Matrix<T,1>::slice(Index n, Index m) const
// m elements starting with a[n]
{
	if (n<0) n=0;
	else if(d1<n) n=d1;    // one beyond the end
	if (m<0 || (d1<n+m)) m=d1-n;
	Row<T,1> z(n,m,*this);
	return z;
}

/** Matrix<T,1>-Matrix<T,1> operations **/
template <class T>
Matrix<T,2> Matrix<T,1>::operator*(const Matrix<T,1>& m1) {
		if((*this).size()==0 || m1.size() == 0) throw Matrix_exception("unitialize matrix in Matrix<T,1> operator*");
		Index d1 = this->dim1();
		Index d2 = m1.dim1();
		if(d1 != d2) throw Matrix_exception("Matrix dimensions don't match in Matrix<T,1> operator*");

		Matrix<T,2> m2(d1,d2);

 		 for(Index i=0;i<d1;++i) {
 		    for(Index k=0;k<d2;++k) {
				m2(i,k) += (*this)(i)*m1(k);
			}
		}
	return m2;
}
template <class T>
Matrix<T,1> Matrix<T,1>::operator+(const Matrix<T,1>& m1) {
		if((*this).size()==0 || m1.size() == 0) throw Matrix_exception("unitialize matrix in Matrix<T,1> operator+");
		Index d1 = this->dim1();
		if(d1 != m1.dim1()) throw Matrix_exception("Matrix dimensions don't match in Matrix<T,1> operator+");
		Index i = 0;
		Matrix<T,1> m2(d1);
		
		for(i=0;i<d1;++i)
				m2(i)=(*this)(i)+m1(i);
		return m2;
}

template <class T>
Matrix<T,1> Matrix<T,1>::operator-(const Matrix<T,1>& m1) {
		if((*this).size()==0 || m1.size() == 0) throw Matrix_exception("unitialize matrix in Matrix<T,1> operator-");
		Index d1 = this->dim1();
		if(d1 != m1.dim1()) throw Matrix_exception("Matrix dimensions don't match in Matrix<T,1> operator-");
		Index i = 0;
		Matrix<T,1> m2(d1);
		
		for(i=0;i<d1;++i)
				m2(i)=(*this)(i)-m1(i);
		return m2;
}

template <class T>
Matrix<T,1> Matrix<T,1>::operator&(const Matrix<T,1>& m1) {
		if((*this).size()==0 || m1.size() == 0) throw Matrix_exception("unitialize matrix in Matrix<T,1> operator&");
		Index d1 = this->dim1();
		if(d1 != m1.dim1()) throw Matrix_exception("Matrix dimensions don't match in Matrix<T,1> operator&");
		Index i = 0;
		bool result = true;
		Matrix<T,1> m2(d1);
		
		for(i=0;i<d1;++i)
				m2(i) = (*this)(i)&m1(i);
		return m2;
}

template <class T>
Matrix<T,1> Matrix<T,1>::operator|(const Matrix<T,1>& m1) {
		if((*this).size()==0 || m1.size() == 0) throw Matrix_exception("unitialize matrix in Matrix<T,1> operator|");
		Index d1 = this->dim1();
		if(d1 != m1.dim1()) throw Matrix_exception("Matrix dimensions don't match in Matrix<T,1> operator|");
		Index i = 0;
		bool result = true;
		Matrix<T,1> m2(d1);
		
		for(i=0;i<d1;++i)
				m2(i) = (*this)(i)|m1(i);
		return m2;
}

template <class T>
Matrix<T,1> Matrix<T,1>::operator^(const Matrix<T,1>& m1) {
		if((*this).size()==0 || m1.size() == 0) throw Matrix_exception("unitialize matrix in Matrix<T,1> operator^");
		Index d1 = this->dim1();
		if(d1 != m1.dim1()) throw Matrix_exception("Matrix dimensions don't match in Matrix<T,1> operator^");
		Index i = 0;
		bool result = true;
		Matrix<T,1> m2(d1);
		
		for(i=0;i<d1;++i)
				m2(i) = (*this)(i)^m1(i);
		return m2;
}

/******************************* Matrix<T,2> **********************************/
template <class T>
Matrix<T,2>& Matrix<T,2>::operator=(const Matrix<T,2>& a)
// copy assignment
{
	if((*this).size()==0) throw Matrix_exception("copying into uninitialized matrix in operator=");
	if(a.size()==0) throw Matrix_exception("copying from uninitialized matrix in operator=");
	if (d1!=a.d1 || d2!=a.d2) throw Matrix_exception("length error in 2D in operator=");
		this->base_assign(a);
	return *this;
}

template <class T>
Matrix<T,2>& Matrix<T,2>::operator=(Matrix<T,2>&& a)
// move assignment
{
	xfer(a);
	(*this).elem=std::move(a.elem);
	(*this).iter=(*this).elem.begin();
	return *this;
}

template <class T>
Row<T,2> Matrix<T,2>::slice(Index n1, Index n2, Index length1, Index length2)
// rows [n:d1)
{
	if (n1<0) n1=0;
	else if(d1<n1) n1=d1;  // one beyond the end
	if (length1<0 || (d1<n1+length1)) length1=d1-n1;
	
	if (n2<0) n2=0;
	else if(d2<n2) n2=d2;   // one beyond the end
	if (length2<=0 || (d2<n2+length2)) length2=d2-n2;

	Row<T,2> z(n1,n2,length1,length2,(*this));
	return z;
}

template <class T>
const Row<T,2> Matrix<T,2>::slice(Index n1, Index n2, Index length1, Index length2) const
// rows [n:d1)
{
	if (n1<0) n1=0;
	else if(d1<n1) n1=d1;  // one beyond the end
	if (length1<0 || (d1<n1+length1)) length1=d1-n1;
	
	if (n2<0) n2=0;
	else if(d2<n2) n2=d2;   // one beyond the end
	if (length2<=0 || (d2<n2+length2)) length2=d2-n2;

	Row<T,2> z(n1,n2,length1,length2,(*this));
	return z;
}


template <class T>
Matrix<T,2> Matrix<T,2>::operator*(const Matrix<T,2>& m1) {
		T temp_A;
		if((*this).size()==0 || m1.size() == 0) throw Matrix_exception("unitialize matrix in Matrix<T,2> operator*");
		Index d1 = this->dim1();
		Index d2 = this->dim2();
		Index d3 = m1.dim2();
		if(d2 != m1.dim1()) throw Matrix_exception("Matrix dimensions don't match in Matrix<T,2> operator*");
		Index block_size_i = 512;
		Index block_size_k = 512;
		Index block_size_j = 2048;
		Index block_i = 0;
		Index block_j = 0;
		Index block_k = 0;
		Matrix<T,2> m2(d1,d3,0);

 		 for(Index i=0;i<d1;i+=block_size_i) {
 		    for(Index k=0;k<d2;k+=block_size_k) {
 		       for(Index j=0;j<d3;j+=block_size_j) {
			   if((i+block_size_i)<d1) {
			        block_i = (i+block_size_i);
			   }
			   else block_i = d1;
                           if((k+block_size_k)<d2) {
                                block_k = (k+block_size_k);
                           }
                           else block_k = d2;
                           if((j+block_size_j)<d3) {
                                block_j = (j+block_size_j);
                           }
                           else block_j = d3;

        		  for(Index ii=i;ii<block_i;ii++) {
        		     for(Index kk=k;kk<block_k;kk++) {
	              		  temp_A = (*this)(ii,kk);
               			 for(Index jj=j;jj<block_j;jj++)  {
							m2(ii,jj) += temp_A*m1(kk,jj);
						}
					}
				}
   		    }
		}
	}
	return m2;
}


template <class T>
Matrix<T,2> Matrix<T,2>::operator+(const Matrix<T,2>& m1) {
		if((*this).size()==0 || m1.size() == 0) throw Matrix_exception("unitialize matrix in Matrix<T,2> operator+");
		Index l1 = this->dim1();
		Index l2 = this->dim2();
		if(l1 != m1.dim1() || l2 != m1.dim2()) throw Matrix_exception("Matrix dimensions don't match in Matrix<T,2> operator+");
		Index i = 0;
		Index j = 0;
		Matrix<T,2> m2(l1,l2);
		
		for(i=0;i<l1;++i)
			for(j=0;j<l2;++j)
				m2(i,j)=(*this)(i,j)+m1(i,j);
		return m2;;
}

template <class T>
Matrix<T,2> Matrix<T,2>::operator-(const Matrix<T,2>& m1) {
		if((*this).size()==0 || m1.size() == 0) throw Matrix_exception("unitialize matrix in Matrix<T,2> operator-");
		Index l1 = this->dim1();
		Index l2 = this->dim2();
		if(l1 != m1.dim1() || l2 != m1.dim2()) throw Matrix_exception("Matrix dimensions don't match in Matrix<T,2> operator-");
		Index i = 0;
		Index j = 0;
		Matrix<T,2> m2(l1,l2);
		
		for(i=0;i<l1;++i)
			for(j=0;j<l2;++j)
				m2(i,j)=(*this)(i,j)-m1(i,j);
		return m2;
}

template <class T>
Matrix<T,2> Matrix<T,2>::operator&(const Matrix<T,2>& m1) {
		if((*this).size()==0 || m1.size() == 0) throw Matrix_exception("unitialize matrix in Matrix<T,2> operator&");
		Index l1 = this->dim1();
		Index l2 = this->dim2();
		if(l1 != m1.dim1() || l2 != m1.dim2()) throw Matrix_exception("Matrix dimensions don't match in Matrix<T,2> operator&");
		Index i = 0;
		Index j = 0;
		Matrix<T,2> m2(l1,l2);
		
		for(i=0;i<l1;++i)
			for(j=0;j<l2;++j)
				m2(i,j)=(*this)(i,j)&m1(i,j);
		return m2;
}

template <class T>
Matrix<T,2> Matrix<T,2>::operator|(const Matrix<T,2>& m1) {
		if((*this).size()==0 || m1.size() == 0) throw Matrix_exception("unitialize matrix in Matrix<T,2> operator|");
		Index l1 = this->dim1();
		Index l2 = this->dim2();
		if(l1 != m1.dim1() || l2 != m1.dim2()) throw Matrix_exception("Matrix dimensions don't match in Matrix<T,2> operator|");
		Index i = 0;
		Index j = 0;
		Matrix<T,2> m2(l1,l2);
		
		for(i=0;i<l1;++i)
			for(j=0;j<l2;++j)
				m2(i,j)=(*this)(i,j)|m1(i,j);
		return m2;
}

template <class T>
Matrix<T,2> Matrix<T,2>::operator^(const Matrix<T,2>& m1) {
		if((*this).size()==0 || m1.size() == 0) throw Matrix_exception("unitialize matrix in Matrix<T,2> operator^");
		Index l1 = this->dim1();
		Index l2 = this->dim2();
		if(l1 != m1.dim1() || l2 != m1.dim2()) throw Matrix_exception("Matrix dimensions don't match in Matrix<T,2> operator^");
		Index i = 0;
		Index j = 0;
		Matrix<T,2> m2(l1,l2);
		
		for(i=0;i<l1;++i)
			for(j=0;j<l2;++j)
				m2(i,j)=(*this)(i,j)^m1(i,j);
		return m2;
}


/******************************* Matrix<T,3> **********************************/
template <class T>
Matrix<T,3>& Matrix<T,3>::operator=(const Matrix<T,3>& a)
// copy assignment
{
	if((*this).size()==0) throw Matrix_exception("copying into uninitialized matrix in Matrix<T,2> operator=");
	if(a.size()==0) throw Matrix_exception("copying from uninitialized matrix in Matrix<T,2> operator=");
	if (d1!=a.d1 || d2!=a.d2 || d3!=a.d3) throw Matrix_exception("length error in 3D in Matrix<T,2> operator=");
		this->base_assign(a);
	return *this;
}

template <class T>
Matrix<T,3>& Matrix<T,3>::operator=(Matrix<T,3>&& a)
// move assignment
{
	xfer(a);
	(*this).elem=std::move(a.elem);
	(*this).iter=(*this).elem.begin();
	return *this;
}

template <class T>
Row<T,3> Matrix<T,3>::slice(Index n1, Index n2, Index n3, Index length1, Index length2, Index length3)
// rows [n:d1)
{
	if (n1<0) n1=0;
	else if(d1<n1) n1=d1;  // one beyond the end
	if (length1<=0 || (d1<n1+length1)) length1=d1-n1;
	
	if (n2<0) n2=0;
	else if(d2<n2) n2=d2;   // one beyond the end
	if (length2<=0 || (d2<n2+length2)) length2=d2-n2;

	if (n3<0) n3=0;
	else if(d3<n3) n3=d3;   // one beyond the end
	if (length3<=0 || (d2<n3+length3)) length3=d3-n3;
	
	Row<T,3> z(n1,n2,n3,length1,length2,length3,(*this));
	return z;
}

template <class T>
const Row<T,3> Matrix<T,3>::slice(Index n1, Index n2, Index n3, Index length1, Index length2, Index length3) const
// rows [n:d1)
{
	if (n1<0) n1=0;
	else if(d1<n1) n1=d1;  // one beyond the end
	if (length1<=0 || (d1<n1+length1)) length1=d1-n1;
	
	if (n2<0) n2=0;
	else if(d2<n2) n2=d2;   // one beyond the end
	if (length2<=0 || (d2<n2+length2)) length2=d2-n2;

	if (n3<0) n3=0;
	else if(d3<n3) n3=d3;   // one beyond the end
	if (length3<=0 || (d2<n3+length3)) length3=d3-n3;
	
	
	Row<T,3> z(n1,n2,n3,length1,length2,length3,(*this));
	return z;
}


template <class T>
Matrix<T,3> Matrix<T,3>::operator+(const Matrix<T,3>& m1) {
		if((*this).size()==0 || m1.size() == 0) throw Matrix_exception("unitialize matrix in Matrix<T,3> operator+");
		Index l1 = this->dim1();
		Index l2 = this->dim2();
		Index l3 = this->dim3();
		if(l1 != m1.dim1() || l2 != m1.dim2() || l3 != m1.dim3()) throw Matrix_exception("Matrix dimensions don't match in Matrix<T,3> operator+");
		Index i = 0;
		Index j = 0;
		Index k = 0;
		Matrix<T,2> m2(l1,l2,l3);
		
		for(i=0;i<l1;++i)
			for(j=0;j<l2;++j)
				for(k=0;k<l3;++k)
					m2(i,j,k)=(*this)(i,j,k)+m1(i,j,k);
		return m2;
}

template <class T>
Matrix<T,3> Matrix<T,3>::operator-(const Matrix<T,3>& m1) {
		if((*this).size()==0 || m1.size() == 0) throw Matrix_exception("unitialize matrix in Matrix<T,3> operator-");
		Index l1 = this->dim1();
		Index l2 = this->dim2();
		Index l3 = this->dim3();
		if(l1 != m1.dim1() || l2 != m1.dim2() || l3 != m1.dim3()) throw Matrix_exception("Matrix dimensions don't match in Matrix<T,3> operator-");
		Index i = 0;
		Index j = 0;
		Index k = 0;
		Matrix<T,2> m2(l1,l2,l3);
		
		for(i=0;i<l1;++i)
			for(j=0;j<l2;++j)
				for(k=0;k<l3;++k)
					m2(i,j,k)=(*this)(i,j,k)-m1(i,j,k);
		return m2;
}
template <class T>
Matrix<T,3> Matrix<T,3>::operator&(const Matrix<T,3>& m1) {
		if((*this).size()==0 || m1.size() == 0) throw Matrix_exception("unitialize matrix in Matrix<T,3> operator&");
		Index l1 = this->dim1();
		Index l2 = this->dim2();
		Index l3 = this->dim3();
		if(l1 != m1.dim1() || l2 != m1.dim2() || l3 != m1.dim3()) throw Matrix_exception("Matrix dimensions don't match in Matrix<T,3> operator&");
		Index i = 0;
		Index j = 0;
		Index k = 0;
		Matrix<T,2> m2(l1,l2,l3);
		
		for(i=0;i<l1;++i)
			for(j=0;j<l2;++j)
				for(k=0;k<l3;++k)
					m2(i,j,k)=(*this)(i,j,k)&m1(i,j,k);
		return m2;
}

template <class T>
Matrix<T,3> Matrix<T,3>::operator|(const Matrix<T,3>& m1) {
		if((*this).size()==0 || m1.size() == 0) throw Matrix_exception("unitialize matrix in Matrix<T,3> operator|");
		Index l1 = this->dim1();
		Index l2 = this->dim2();
		Index l3 = this->dim3();
		if(l1 != m1.dim1() || l2 != m1.dim2() || l3 != m1.dim3()) throw Matrix_exception("Matrix dimensions don't match in Matrix<T,3> operator|");
		Index i = 0;
		Index j = 0;
		Index k = 0;
		Matrix<T,2> m2(l1,l2,l3);
		
		for(i=0;i<l1;++i)
			for(j=0;j<l2;++j)
				for(k=0;k<l3;++k)
					m2(i,j,k)=(*this)(i,j,k)|m1(i,j,k);
		return m2;
}

template <class T>
Matrix<T,3> Matrix<T,3>::operator^(const Matrix<T,3>& m1) {
		if((*this).size()==0 || m1.size() == 0) throw Matrix_exception("unitialize matrix in Matrix<T,3> operator^");
		Index l1 = this->dim1();
		Index l2 = this->dim2();
		Index l3 = this->dim3();
		if(l1 != m1.dim1() || l2 != m1.dim2() || l3 != m1.dim3()) throw Matrix_exception("Matrix dimensions don't match in Matrix<T,3> operator^");
		Index i = 0;
		Index j = 0;
		Index k = 0;
		Matrix<T,2> m2(l1,l2,l3);
		
		for(i=0;i<l1;++i)
			for(j=0;j<l2;++j)
				for(k=0;k<l3;++k)
					m2(i,j,k)=(*this)(i,j,k)^m1(i,j,k);
		return m2;
}

//-----------------------------------------------------------------------------
}
#endif
